PerturbedEquilibrium - FEATURE - Decompose the plasma energy of free-boundary eigenmodes on the (psi,theta) grid - #465

Consolidated review from three independent passes (structural, performance, code clarity), each claim checked against the source. Ordered by what I'd want settled before merge. The physics itself was not reviewed here — see the end.

1. Blocking: the stage aborts the run instead of skipping

run_energy_decomposition (GeneralizedPerturbedEquilibrium.jl, new hunk) gates only require_solution and require(:free_boundary), which warn and return. But decompose_energy throws ArgumentError on nlow != nhigh and mats.kinetic !== nothing (EnergyDecomposition.jl:698–699). A legal multi-n or kinetic deck that carries [EnergyDecomposition] therefore kills the whole run after FFS — before PE, KineticForces or SLAYER — contradicting the stage docstring ("warned and skipped, like every optional stage") and docs/src/energy_decomposition.md. Copilot flagged the same gap. ~10 lines: warn-and-return for both conditions in the driver, keep the throws for scripted use.

2. Decide before merge: where the output lives

• The group ForceFreeStates/EnergyDecomposition/ bakes "free-boundary eigenmodes only" into the schema, while the producer is a PE-module stage. Free-boundary enters the code through just wt[:,k] (:342) and ep[k] (:814); a decompose_energy(ffs, edge::Vector; reference=nothing) entry would serve fixed-boundary modes and forced responses at ~half a day, and PerturbedEquilibrium/… would then be the natural home. HDF5 renames are clean breaks here with no shim, so the placement is worth deciding now even if the generalization waits.
• eigenmode_index = copy(eigenmodes) (:719) echoes a control value into an output group (hdf5-conventions.md: never echo control-struct values outside Input/). It is defensible as the coordinate of the eigenmode dimension — in which case attach it as a dimension scale; only psi is attached today (:926–927).

3. Structural: the geometry block shadows Equilibrium

• surface_geometry! (:160–166) reproduces Equilibrium.flux_surface_metric (FluxSurfaceMetrics.jl:16–23) line for line — the function whose docstring declares it "the single source of truth for the |∇ψ| computation reused across the FFS and PE modules" — then extends it with w21/w22/w31/w32 and the tangent basis. Extend flux_surface_metric (or add a sibling in Equilibrium) and call it, rather than carrying a second implementation in PE.
• The same block re-interpolates geometry with 13 bicubic evaluations per point at the equilibrium knots, where rzphi_*.nodal_derivs.partials[1:3, i, j] already holds f, ∂f/∂ψ, ∂f/∂θ and make_metric (Fourfit.jl:193–197) and GeometryProfiles.jl:44 read them directly. Only the appended edge node needs the interpolant. The feature is cheap overall (order 0.1–1 s), so this is about consistency with the existing pattern more than speed.
• metric_bsq_interpolant (:213–230) is a private, persistent B² that competes with equil.eqfun_B and is rebuilt on every call (:717). Whether metric-consistent B² is the right physics is for the physics review; either way its owner should be Equilibrium or make_metric, so a future consumer (LocalStability shear, KineticForces) does not grow a third B².
• eigenmode_modes (:361–367) re-encodes −A⁻¹(Bξ′ + Cξ) by hand; compute_node_xi_s! (ForceFreeStates/EulerLagrange.jl:1397, used at Result.jl:174) already does exactly that.
• surface_geometry! and curvature_kernel! depend only on (equil, psi) but sit inside the eigenmode loop (:755–770): all geometry is computed twice for eigenmodes = [1, 2]. Hoist per surface.

4. Robustness

• surface_geometry!(…; bsq=nothing) leaves geom.bsq zero, and curvature_kernel! divides by it (:285). Make bsq required.
• similar(x, 0) (:921–922) collapses disabled 2-D/3-D outputs to rank-1 while the metadata table declares (psi, eigenmode) / (psi, theta, eigenmode) — the writer cannot attach the declared dims (Copilot). The zero-extent sentinel has precedent, but only if the rank is kept.
• Single-n is load-bearing, not just declared: N = ffs.numpert_total (:349) is written into xi_s[i, :] of width mpert (:367), a DimensionMismatch the guard at :698 is what prevents.
• Constructors hard-code field counts (for _ in 1:16, 1:6, 1:13, 1:6, 1:8 at :131/248/400/451/539); fieldcount(T) stops them rotting when a field is added.

5. Cleanup

• Code comments are clean — no run exceeds two lines. The narration is in docstrings: :208–211 and :409–412 carry session history (keep the one-sentence reason; drop the rest), and test:176, :208 record observed percentages that belong in the docs, not a comment.
• docs/src/energy_decomposition.md:64 says "better than 1 %"; the tests record 1.4 %, 1.2 % and 3.3 % at the mpsi = 16 fixture (test:176, :208). "≲1–3 % at mpsi = 16, converging with resolution" matches the evidence.
• test/runtests_energy_decomposition.jl (292 lines) would be ~200 with the same coverage: _ed_metadata_violations (:216–231) duplicates test/h5_metadata_check.jl — include it; :32–50 tests the reg_spot keyword, not this kernel; :70–73 asserts kwdef defaults; :86–90 re-derive K1/K2/K3 from :293–295 (tautological); :99–128 hand-plumbs the curl residual that :173 already pins; decompose_energy runs eight times — hoist the shared results. :287–291 writes into the shared fixture's gpec.h5 and only works because it is the last testset — write to a temp copy.
• Five structs describe their fields in prose rather than a ## Fields block (:104, :232, :376, :435, :516); the unnamed 1e-10 at :315; μ₀j^θ and μ₀j^ζ computed three times (:286, :475, :572); the 13 hand-named accumulators plus three mirrored copy blocks (:737–749, :796–837) could be driven from one name tuple, as _allocate_result! already half-does.

6. Before merge

• The branch is ~250 commits behind and develop has since rewritten FieldReconstruction.jl (225 lines) and the driver (+455); your 6-line reg_spot change and 47-line hook land on rewritten code. DIRTY here is semantic, not textual (zero conflict markers) — worth re-reading both hunks after the merge-forward.
• The regression report was a laptop run (arm64, 4 threads) against a develop that has moved; please re-run on the shared machine after merging forward and refresh the stamp.
• This review deliberately did not cover the physics (metric-B² vs eqfun_B, the pointwise covariant lowering your docstring at :408–412 justifies, the ξ_s recomputation at knots). Those deserve a physics pass and I'd ask for one before approval.

Credit where due: Copilot's two inline comments — the gating gap and the similar(x, 0) rank collapse — were both real and are folded in above.

🤖 Generated with Claude Code

https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU

Physics pass, to close the gap my earlier review left open. Checked term by term against the Fortran counterpart and the theory note on origin/recon of the Fortran repository (gpout_recon/gpout_recon2 in gpout.f, gpeq_c/gpeq_K/gpeq_shear/gpeq_curvature/gpeq_cveri/gpeq_firstform in gpeq.f, docs/tex/recon/main.tex) — not on Fortran develop, which is worth stating in the docs since "fidelity to Fortran" means fidelity to that branch. Every claim marked verified was read or derived; the Fortran line references were checked by me as well.

Verdict: correct with caveats. No missing terms, no placeholders; every kernel either matches origin/recon line for line or is a defensible, more accurate deviation. One thing should be settled before merge — see §6 — because it is the untested explanation for the 1.4–3.3 % mismatch.

1. Both δW forms — verified

• Standard form (:565–585): derived from −½∫ξ*·F(ξ) with the full ξ (the ξ∥ pieces of the J×Q and ∇·ξ terms cancel, so using the full Jξ^θ, Jξ^ζ is right). det_term at :572–575 is J·det[j, b, ξ*] with j^ψ = 0; ∇·ξ at :577 is (1/J)[∂ψ(Jξ^ψ) + ∂θ(Jξ^θ) − 2πin Jξ^ζ] in the exp(2πi(mθ−nζ)) convention. Term for term identical to gpeq_firstform. One divergence: ∂θ(Jξ^θ) is spectral here (:415–418) and a periodic spline there — fine, but annotate it as a Fortran deviation.
• Effective field b_eff = b + ξ_n μ₀ j×n̂ (:455–480): all six component formulas verified by hand and identical to gpeq_c (cwt/cwz/cvp/cvt/cvz). μ₀j^θ = −F′/J and μ₀j^ζ = −(μ₀p)′/χ′ − qF′/J follow from ∇×B and (j×B)·∇ψ = p′, given F_deriv = d(2πF)/dψ_N (inferred from the Fortran sq%f1(1) usage; please confirm).
• K = K1 + K2 + K3 (:293–295): identical to gpeq_K and main.tex; units consistent (K|ξ_n|² is Pa); K3 = 2p′κ_ψ checks against the interchange term. Caveat: K1, K2 and the completing-the-square identity were verified against main.tex and gpeq_K only — Bernstein 1983, Chance 1992 and Boozer 2004 are not in docs/resources/, so the equation numbers this PR adds to references.md/citations.md are unverified. Please either add the PDFs or cite the theory note.
• ½, Re, conjugation, normalization: Re[(Jc^i)* c_i]/(μ₀J) is Hermitian-correct; the ½ matches Fortran dw_raw = 0.5*total; the 2μ₀/ψ₀² scale (:602) with unit-power-normalized wt[:,k] as edge BC (:342) is the same convention as Fortran's dcon_fac after idcon.f undoes its rescaling; ep = diag(wt†·wp·wt) with wp = (u2/u1)/ψ₀² matches free.f. No vacuum or surface term belongs in the comparison — ep is the plasma boundary term Ξ†(FΞ′+KΞ), i.e. the volume δW_p; ev is separate. Verified in concept; the DCON ½/μ₀ scale was not re-derived from Glasser 2016, but the 1e-3 agreement rules out a missing O(1) factor.

2. Metric-B² vs eqfun_B — your choice is right; say so at the call site

Fortran gpeq_K and gpeq_curvature take B² from eqfun (bsq_val = eqfun%f(1)**2, with fx/fy for its derivatives). make_metric builds only g_ij/J, J, J_ψ — no eqfun — so the Euler–Lagrange operator and the b reconstruction carry the metric B², and the C∥ = K2|ξ_n|² cancellation holds only with that B². Using it is operator-consistent and an improvement over Fortran. The ~1 % is the mismatch between the equilibrium's direct |B| and the inverse-map metric at mpsi = 16 — coordinate-transformation accuracy converging with equilibrium resolution, not hidden physics. It needs the project's one-line deviation note at :285, :291, :294, not only in the docstring.

3. Pointwise covariant lowering — the claim holds, with two corrections

gpeq_c back-transforms band-limited modes before squaring, so ∮|C|² and ∮K2|ξ_n|² each lose out-of-band energy that does not cancel between them; pointwise products on the mtheta grid are Parseval-exact (aliasing aside). So the Julia effective-field numbers will differ from Fortran recon reference values at the few-% level, and Julia is the more accurate. Two corrections to the docstring at :408–412: "resolution-independent" should read "ψ-resolution-independent, mpert-dependent" (the deficit shrinks with the band); and this is a hybrid — Jξ^ψ is formed pointwise (:416, :426) but Jξ^θ, Jξ^ζ still come from the truncated jmat convolution in compute_contra_displacements — the docstring implies all-pointwise.

4. ξ_s at the knots — verified, no physical content in knot-vs-node

−A⁻¹(Bξ′ + Cξ) (:361–367) matches Fortran gpeq_sol and Julia's own compute_node_xi_s! (EulerLagrange.jl:1397, same raw A; Cholesky vs LDLᵀ immaterial). It is the algebraic Euler–Lagrange relation for Ξ_s, so evaluating at knots rather than ODE nodes changes only the interpolation accuracy of (ξ, ξ′), nothing physical. That argues for calling the existing helper rather than re-encoding it (as my earlier review noted).

5. The curl residual — right check, state its scope

(:495–515) is identical to gpeq_cveri and is the Ξ_s equation itself (∇ψ·∇×C = 0 ⇔ δj·∇ψ = 0 closure), so at the knots it holds by construction. It validates the ξ_s algebra, the g22/g23/g33 lowering and the b_θ, b_ζ assembly against A, B, C — not the ψ-ODE and not b_ψ or g11/g12/g13/|∇ψ|². A roundoff-level value says nothing about the eigenmode's correctness. Worth one sentence in the docs so it isn't read as a global check.

6. Convergence — the one thing to settle before merge

Expected: ~h⁴ for the cubic-spline ψ-quadrature on a smooth integrand, spectral in θ. Observed from the PR's own numbers (mpsi 16 → 128): effective form ~h^0.9, standard form ~h^1.8. That is not quadrature-order behaviour.

The leading explanation is a mis-documented Fortran deviation. gpout.f says verbatim, just above the energy check: "Perform the energy check on the DCON solution grid. Do not use rzphi%xs for this comparison; it is a coarse geometry grid." — and integrates on psifac. decomposition_grid (:309–316) integrates on the rzphi_xs equilibrium knots, and the docstring at :305 says gpout_recon "loops over the same knots" — true of its diagnostics loop, not of its energy integral. With the 1/(m − nq) structure in ξ_s and Jξ^θ, the coarse-knot integral is the natural suspect for both the size of the mismatch and the sub-quadrature rate; the slower effective-form rate is consistent with K depending on spline derivatives (bsq_psi, dratio) that lose an order. A confound: ep itself moves with mpsi.

Discriminating test, cheap to run: fix mpsi, integrate on ODE nodes ∪ knots instead of knots alone, and report how the relative error changes. If it collapses, fix the grid and the docstring; a separate mpert scan would then bound the band-truncation floor at mpsi = 128.

7. Fortran correspondence and the annotations owed

gpout_recon/gpeq_c/gpeq_K/gpeq_shear/gpeq_curvature/gpeq_cveri/gpeq_epf/gpeq_dst/gpeq_firstform ↔ decompose_energy/effective_field!/curvature_kernel!/curl_residual!/surface_densities!. Term-level divergences, each needing the project's one-line deviation comment at the call site: (a) metric B² vs eqfun B²; (b) pointwise vs band-truncated C_i, ξ_n; (c) spectral vs spline ∂θ; (d) radial grid knots vs psifac; (e) reg_spot = 0 vs Fortran's regularized xmp1/xmt path in gpeq_firstform (inferred from variable names, not traced). Also: the long_name on current_coupling and pressure_compression (:940–941) says "imaginary part is roundoff" — true of the volume totals, not of the per-surface profiles; reword.

Not checked: the Bernstein/Chance/Boozer equation numbers; the DCON ½/μ₀ convention from Glasser 2016 directly; the edge amplitude in sum_eigenmode_contributions; gpeq_firstform's regularization path; gpout_recon2's radial grid; aliasing error of the mtheta Riemann sum.

🤖 Generated with Claude Code

https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU